ESA18_Workshop_AusCover

NORMALIZED DIFFERENCE VEGETATION INDEX (NDVI) CHANGE ANALYSIS

In this section we will explore the change in NDVI in an area west of Tara where Coal Seam Gas development has been happening in recent years. NDVI can be used to estimate the density of green vegetation on a patch of land. It is derived from remote sensing data, typically a space platform. Chlorophyll, the pigment in plant leaves, strongly absorbs red and blue light. Leaves, on the other hand, strongly reflect near-infrared light (NIR) and green light (this is why vegetation looks green to our eyes). NDVI takes advantage of this difference in behaviour for different light wavelengths by green vegetation, combining measurements at different wavelengths in this formula: NDVI = (NIR-Red)/(NIR+Red). NDVI ranges from -1 to 1. With different values indicating different situations:

Values close to ‘+1’ (0.8-0.9): Very high density of green leaves.
Low values: Low vegetation.
Very low values (e.g. 0.2): Nearly bare soils.
Values close to ‘0’: Practically no vegetation (e.g. densely populated urban areas).
Very low positive to slightly negative values: Free standing water (e.g. oceans, seas, lakes, rivers), which have rather low reflectance in both spectral bands.
Values close to ‘-1’: likely to represent atmospheric disturbance or some specific materials.

Generally, to investigate NDVI change over time an atmospheric correction needs to be performed. The data that we use in this workshop has already been corrected to reduce contamination by cloud and other problems.

Preparation: Getting ready

We start by loading the required libraries, setting up a working directory, and cleaning up memory.

# Load  Libraries
# ===============

library(dplyr)
library(reshape2)
library(stringr)

library(RColorBrewer)
library(ggplot2)
library(gridExtra)
library(RStoolbox)

library(sp)

library(rgdal)
library(raster)
library(rasterVis)

##library(RgoogleMaps)
library(ggmap)

#library(maps)
#library(mapdata)
#library(maptools)


# Set Data Path (TERN AusCover data repository)
# =============================================
data.path = "http://qld.auscover.org.au/public/data/landsat/surface_reflectance/aus"


# Optional Steps (remove comments to run them)
# ==============

# Setting up my Current Working Directory
# ---------------------------------------
#getwd()
#my.CWDir = "C:/Users/uqbblanc/Documents/TERN/CWDir"
#setwd(my.CWDir)
#getwd()

# Clean up Memory
# ---------------
#rm(list=ls())
#list.files()

Explore an Area where NDVI Change might have occurred

We import and display a Stamen map to explore the area that we will be working with (West of Tara). The first attempt, using the extent of this area, leads to too detailed map. We add a plot where we zoom out by 2 degrees and now we can identify a nearby location (Dalby).

# Define the Extent of the Area to Explore
# ----------------------------------------
WofTara.extent = extent(150.3727, 150.4963, -27.0767, -26.9819)
WofTara.extent

## class       : Extent 
## xmin        : 150.3727 
## xmax        : 150.4963 
## ymin        : -27.0767 
## ymax        : -26.9819

# Get Map with for West of Tara Extent 
# ------------------------------------
WofTara.StamenMap = get_stamenmap(WofTara.extent[c(1,3,2,4)], maptype="terrain")

# Create Plot of the Map
# ----------------------
WofTara.Map.P1 = 
ggmap(WofTara.StamenMap) + labs(title= "Study Area (W of Tara): Too detailed", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"), 
      axis.title = element_text(size=11, face="bold"), axis.text=element_text(size=9) )
#WofTara.Map.P1  # Plot individually for a bigger image

# Zoom out by 2 degrees & Create New Plot
# ---------------------------------------
 # Download Map with new wider extent
WofTaraPlus2.StamenMap = get_stamenmap((WofTara.extent+2)[c(1,3,2,4)], maptype="terrain")
 # Create Plot for new wider map
WofTaraPlus2.Map.P2 = 
ggmap(WofTaraPlus2.StamenMap) + labs(title= "Study Area (W of Tara) + 2 Degrees", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"), 
      axis.title = element_text(size=11, face="bold"), axis.text=element_text(size=9) )
 # Create Polygon with Original Extent and Add it to the Plot
WofTara.SP = as(WofTara.extent, "SpatialPolygons")
WofTara.SP  # Has no CRS

## class       : SpatialPolygons 
## features    : 1 
## extent      : 150.3727, 150.4963, -27.0767, -26.9819  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

proj4string(WofTara.SP) = CRS("+init=epsg:4326") # Add CRS
WofTara.SP

## class       : SpatialPolygons 
## features    : 1 
## extent      : 150.3727, 150.4963, -27.0767, -26.9819  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

 # Add Polygon to the Plot
WofTaraPlus2.Map.P2 = WofTaraPlus2.Map.P2 + geom_polygon(data=WofTara.SP, aes(x=long, y=lat), color="red", alpha=0)
#WofTaraPlus2.Map.P2  # Plot individually for a bigger image

# Show Both Plots 
# ----------------
grid.arrange(WofTara.Map.P1, WofTaraPlus2.Map.P2, nrow=1)

Create a function to Download, Read, Subset, and Combine the Landsat Surface Reflectance Data.

This function includes many of the steps in the INTRODUCTORY STEPS section. Specifically the function perform these tasks:

Download the data file
Load the data file required bands
Subset the RasterLayer objects
Replaced values < 0 with the Missing Data code (Not really required in this example)
Create a multi-layered RasterBrick object
Return the multi-layered RasterBrick object

The function takes the following arguments:

data.fpn: Data file path and name
dl.fn: Name for the downloaded file
rB.extent: RasterBrick extent
rB.crs: RasterBrick coordinates reference system

There are a number of print commands within the function ‘commented out’. These commands are not required for the correct functioning of the function. However, they have been included because they might make it easier to understand what the function is doing. To investigate further what the function is doing at each step simply remove the comment (#) symbol in front of the relevant print command.

get.B3B4B5.rB.f = function(data.fpn, dl.fn, rB.extent, rB.crs){

    # Download the data file
    # ======================    
    download.file(url=data.fpn, destfile=dl.fn, method='auto')
    #print(list.files())

    
    # Load the data file required bands
    # =================================
    LSSR.Band3.rL = raster(dl.fn, band=3)
    #print(LSSR.Band3.rL)
    LSSR.Band4.rL = raster(dl.fn, band=4)
    #print(LSSR.Band4.rL)
    LSSR.Band5.rL = raster(dl.fn, band=5)
    #print(LSSR.Band5.rL)

    
    # Subset Rasters
    # ==============    

    # Project Extent to a CRS=EPSG:3577
    # ---------------------------------
    #print(rB.extent)
    # Create a Polygon from Extent, Reproject the Poloygon and use its extension
    rB.extent.GeoCoord.SP = as(rB.extent, "SpatialPolygons")
    #print(rB.extent.GeoCoord.SP)  # Has no CRS
    proj4string(rB.extent.GeoCoord.SP) = CRS("+init=epsg:4326")
    #print(rB.extent.GeoCoord.SP)
    rB.extent.EPSG3577.SP = spTransform(rB.extent.GeoCoord.SP, CRS("+init=epsg:3577"))
    #print(rB.extent.EPSG3577.SP)
    rB.extent.EPSG3577 = extent(rB.extent.EPSG3577.SP)
    #print(rB.extent.EPSG3577)

    # Crop (i.e. subset) the raster layers to our area of interest
    # ---------------------------------
    LSSR.Subset.Band3.rL = crop(LSSR.Band3.rL, rB.extent.EPSG3577)
    #print(LSSR.Subset.Band3.rL)
    LSSR.Subset.Band4.rL = crop(LSSR.Band4.rL, rB.extent.EPSG3577)
    #print(LSSR.Subset.Band4.rL)
    LSSR.Subset.Band5.rL = crop(LSSR.Band5.rL, rB.extent.EPSG3577)
    #print(LSSR.Subset.Band5.rL)
    
    # Replace values < 0 with NAs
    # ===========================
     # Band 3: Red
     # ------------
    LSSR.Subset.Band3.rL[LSSR.Subset.Band3.rL < 0] = NA
    #print(LSSR.Subset.Band3.rL)
    #print(freq(is.na(LSSR.Subset.Band3.rL)))
    #print(freq(is.na(LSSR.Subset.Band3.rL))/ncell(LSSR.Subset.Band3.rL))
     # Band 4: NIR
     # ------------
    LSSR.Subset.Band4.rL[LSSR.Subset.Band4.rL < 0] = NA
    #print(LSSR.Subset.Band4.rL)
    #print(freq(is.na(LSSR.Subset.Band4.rL)))
    #print(freq(is.na(LSSR.Subset.Band4.rL))/ncell(LSSR.Subset.Band4.rL))
     # Band 5: SWIR
     # ------------
    LSSR.Subset.Band5.rL[LSSR.Subset.Band5.rL < 0] = NA
    #print(LSSR.Subset.Band5.rL)
    #print(freq(is.na(LSSR.Subset.Band5.rL)))
    #print(freq(is.na(LSSR.Subset.Band5.rL))/ncell(LSSR.Subset.Band5.rL))

    
    # Create a multi-layered RasterBrick object
    # =========================================
    LSSR.Subset.Bands3to5.rB = brick(LSSR.Subset.Band3.rL, LSSR.Subset.Band4.rL, LSSR.Subset.Band5.rL)
    names(LSSR.Subset.Bands3to5.rB) = c("Band3.Red", "Band4.NIR", "Band5.SWIR")
    #print(LSSR.Subset.Bands3to5.rB)


    # Return RasterBrick with the Required Bands Subset to the desired Extent
    # =======================================================================
    return(LSSR.Subset.Bands3to5.rB)

} # get.B3B4B5.rB.f = function(data.fpn, dl.fn, rB.extent, rB.crs){

Obtain Landsat Surface Reflectance (LSSR) Data for the required Extent and Period using the new function (`get.B3B4B5.rB.f`) & Plot them

Use our new function to get the Landsat Surface Reflectance datasets for Autumn 2007 and Autumn 2017, then subset the required bands and finally create and return a RasterBrick with the resulting RasterLayers. The RasterBricks for these periods are then plotted together.

getwd()

## [1] "C:/Users/uqbblanc/Documents/TERN/CWDir/FinalVersions_ToSend/SuperFinal/AusCover"

# Data Path and Data (Landsat Surface Refectance, LSSR) Files Names
data.path = "http://qld.auscover.org.au/public/data/landsat/surface_reflectance/aus"
data.2007.fn = "lztmre_aus_m200703200705_dbia2.vrt"
data.2017.fn = "l8olre_aus_m201703201705_dbia2.vrt"


# Create RasterBricks with LSSR Data for the required Extent and Period
# =====================================================================

# Autumn 2007
# -----------

# Create RasterBrick with data
wofTara.2007.LSSR.rB = get.B3B4B5.rB.f( data.fpn=paste(data.path, data.2007.fn, sep="/"), dl.fn="LSSR_2007.vrt", rB.extent=WofTara.extent, rB.crs=CRS("+init=epsg:4326") )

wofTara.2007.LSSR.rB

## class       : RasterBrick 
## dimensions  : 408, 451, 184008, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : Band3.Red, Band4.NIR, Band5.SWIR 
## min values  :       270,       938,        432 
## max values  :      3404,      4203,       5016

# Create Plot of RasterBrick
wofTara.2007.LSSR.p = ggRGB(wofTara.2007.LSSR.rB, r=3, g=2, b=1) + 
labs(title= "LSSR: Autumn 2007", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=20, face="bold"), 
      axis.title = element_text(size=12, face="bold"), axis.text=element_text(size=9) )

# Autumn 2017
# -----------

# Create RasterBrick with data
wofTara.2017.LSSR.rB = get.B3B4B5.rB.f( data.fpn=paste(data.path, data.2017.fn, sep="/"), dl.fn="LSSR_2017.vrt", 
                                        rB.extent=WofTara.extent, rB.crs=CRS("+init=epsg:4326") )

wofTara.2007.LSSR.rB

## class       : RasterBrick 
## dimensions  : 408, 451, 184008, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : Band3.Red, Band4.NIR, Band5.SWIR 
## min values  :       270,       938,        432 
## max values  :      3404,      4203,       5016

# Create Plot of RasterBrick
wofTara.2017.LSSR.p = ggRGB(wofTara.2017.LSSR.rB, r=3, g=2, b=1) + 
labs(title= "LSSR: Autumn 2017", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=20, face="bold"), 
      axis.title = element_text(size=12, face="bold"), axis.text=element_text(size=9) )

      
# Display both RasterBrick Plots in one plot
# ==========================================
grid.arrange(wofTara.2007.LSSR.p, wofTara.2017.LSSR.p, nrow=1)

Compute NDVI from the extracted Subsets

In R raster calculations can be performed using:

'Raster Algebra' (also known s 'Raster Maths'): Directly combining raster objects in a formula with the required calculations. We can use this formula directly, or place it in a function and call this function. Typically used for small rasters and simple calculations (not recommended otherwise).
'Higher Level Functions': First we place the formula with the required raster calculations in a function and then we use a Higher Level Function to invoke this function on relevant raster objects. This is more efficient (faster) option, which should be used for larger rasters and/or complex formulas. The Higher Level Function used depends on the number of raster object involved in the calculations:
- For a Single Raster Object: calc and overlay are most commonly used.
- For Multiple Raster Objects: overlay is typically used.

See TERN’s DSDP Tutorial ‘Using Raster Data in R’ for details.

To compute NDVI values we first create a formula with the required maths (for details on the NDVI see text at the beginning of this section (NORMALIZED DIFFERENCE VEGETATION INDEX (NDVI) CHANGE ANALYSIS). We actually create two formulas, the first one passing the RasterBrick to be used with calc and the second one passing individual layers to be used with overlay. Then we calculate this index in different ways:

using raster algebra (i.e. without calling any function)
calling the a function we wrote to compute NDVI
using Higher level Functions that call the functions we wrote to calculate the NDVI: (3a) first using calc and (3b) then using overlay.

Sometimes calc and overlay don’t work; we show a trick (generating an empty raster with the right features and initialising it before assigning the values we need for our work to it) to help it work for us. The we compare the results of these 4 approaches to prove that the all produce the same results.

NOTE that in these datasets the value 32767 is used as a no-data code, so we will replace this value by NA

#=========================================================================================
# Create functions to Calculate NDVI
#=========================================================================================

# Function 1: Passing the RasterBrick
# -----------------------------------

calc.NDVI.f1 = function(rb) {

    NDVI.rL = (rb[["Band4.NIR"]] - rb[["Band3.Red"]]) / (rb[["Band4.NIR"]] + rb[["Band3.Red"]])
    return(NDVI.rL)

} # calc.NDVI.f1 = function(rb) {


# Function 2: Passing both RaterLayers
# -----------------------------------

calc.NDVI.f2 = function(rl.B4, rl.B3) {

    NDVI.rL = (rl.B4 - rl.B3) / (rl.B4 + rl.B3)
    return(NDVI.rL)

} # calc.NDVI.f2 = function(rl.B4, rl.B3) {


#=========================================================================================
# (1) Calculate NDVIs using Raster Algebra (= Raster Maths)
#=========================================================================================

# NDVI for 2007
# -------------
wofTara.2007.NDVI.rL.1 = (wofTara.2007.LSSR.rB[["Band4.NIR"]] - wofTara.2007.LSSR.rB[["Band3.Red"]]) / 
                         (wofTara.2007.LSSR.rB[["Band4.NIR"]] + wofTara.2007.LSSR.rB[["Band3.Red"]])
wofTara.2007.NDVI.rL.1

## class       : RasterLayer 
## dimensions  : 408, 451, 184008  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -0.08312343, 0.7602862  (min, max)

# NDVI for 2017
# -------------        
wofTara.2017.NDVI.rL.1 = (wofTara.2017.LSSR.rB[["Band4.NIR"]] - wofTara.2017.LSSR.rB[["Band3.Red"]]) / 
                         (wofTara.2017.LSSR.rB[["Band4.NIR"]] + wofTara.2017.LSSR.rB[["Band3.Red"]])
wofTara.2017.NDVI.rL.1

## class       : RasterLayer 
## dimensions  : 408, 451, 184008  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -0.6510903, 0.8544872  (min, max)

(wofTara.2017.LSSR.rB[["Band3.Red"]] - wofTara.2017.LSSR.rB[["Band4.NIR"]]) / 
                         (wofTara.2017.LSSR.rB[["Band4.NIR"]] + wofTara.2017.LSSR.rB[["Band3.Red"]])

## class       : RasterLayer 
## dimensions  : 408, 451, 184008  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -0.8544872, 0.6510903  (min, max)

#=========================================================================================
# (2) Directly calling one of the functions we created above to compute NDVI
#=========================================================================================
# NDVI for 2007
# -------------
wofTara.2007.NDVI.rL.2 = calc.NDVI.f1(rb = wofTara.2007.LSSR.rB)
wofTara.2007.NDVI.rL.2

## class       : RasterLayer 
## dimensions  : 408, 451, 184008  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -0.08312343, 0.7602862  (min, max)

# NDVI for 2017
# -------------        
wofTara.2017.NDVI.rL.2 = calc.NDVI.f1(rb = wofTara.2017.LSSR.rB)
wofTara.2017.NDVI.rL.2

## class       : RasterLayer 
## dimensions  : 408, 451, 184008  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -0.6510903, 0.8544872  (min, max)

#=========================================================================================
# (3) Using Higher Level Functions
#=========================================================================================

#-----------------------------------------------------------------------------------------
# (3a) Using  'calc' & our 1st function to calculate NDVI
#-----------------------------------------------------------------------------------------
# Section before 'TRICK' commented out because it doesnt' work. Can give it a try
# by removing the comments.

## NDVI for 2007
## .............
#wofTara.2007.NDVI.rL.3a1 = calc(rb=wofTara.2007.LSSR.rB, fun=calc.NDVI.f1)
#wofTara.2007.NDVI.rL.3a1

# NDVI for 2017
# .............
#wofTara.2017.NDVI.rL.3a1 = calc(rb=wofTara.2017.LSSR.rB, fun=calc.NDVI.f1)
#wofTara.2017.NDVI.rL.3a1


# TRICK: Create & Initialise the Raster container and then copy the RasterLayers on it
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# All years have the same dimensions
rL = raster( ncol=ncol(wofTara.2007.LSSR.rB[["Band4.NIR"]]), 
              nrow=nrow(wofTara.2007.LSSR.rB[["Band4.NIR"]]) )
rL.B4.NIR = init(rL, fun=runif)
rL.B3.Red = init(rL, fun=runif)     

# NDVI for 2007
# .............
rL.B4.NIR = wofTara.2007.LSSR.rB[["Band4.NIR"]]
rL.B3.Red = wofTara.2007.LSSR.rB[["Band3.Red"]]
rB = brick(rL.B4.NIR, rL.B3.Red)
names(rB) = c("Band4.NIR", "Band3.Red")
wofTara.2007.NDVI.rL.3a2 = calc(rB, fun=calc.NDVI.f1)
wofTara.2007.NDVI.rL.3a2

## class       : RasterLayer 
## dimensions  : 408, 451, 184008  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -0.08312343, 0.7602862  (min, max)

# NDVI for 2017
# .............
rL.B4.NIR = wofTara.2017.LSSR.rB[["Band4.NIR"]]
rL.B3.Red = wofTara.2017.LSSR.rB[["Band3.Red"]]
rB = brick(rL.B4.NIR, rL.B3.Red)
names(rB) = c("Band4.NIR", "Band3.Red")
wofTara.2017.NDVI.rL.3a2 = calc(rB, fun=calc.NDVI.f1)
wofTara.2017.NDVI.rL.3a2

## class       : RasterLayer 
## dimensions  : 408, 451, 184008  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -0.6510903, 0.8544872  (min, max)

#-----------------------------------------------------------------------------------------
# (3b) Using  'overlay' & our 2nd function to calculate NDVI
#-----------------------------------------------------------------------------------------
# Section before 'TRICK' commented out because it doesnt' work. Can give it a try
# by removing the comments.

## NDVI for 2007
## .............
#wofTara.2007.NDVI.rL.3b1 = overlay( rl.B4=wofTara.2007.LSSR.rB[["Band4.NIR"]],
#                                    rl.B3=wofTara.2007.LSSR.rB[["Band3.Red"]],
#                                    fun=calc.NDVI.f2)
#wofTara.2007.NDVI.rL.3b1

## NDVI for 2017
## .............
#wofTara.2017.NDVI.rL.3b1 = overlay( rl.B4=wofTara.2017.LSSR.rB[["Band4.NIR"]],
#                                    rl.B3=wofTara.2017.LSSR.rB[["Band3.Red"]],
#                                    fun=calc.NDVI.f2)
#wofTara.2017.NDVI.rL.3b1


# TRICK: Create & Initialise the Raster container and then copy the RasterLayers on it
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# All years have the same dimensions
rL = raster( ncol=ncol(wofTara.2007.LSSR.rB[["Band4.NIR"]]), 
              nrow=nrow(wofTara.2007.LSSR.rB[["Band4.NIR"]]) )
rL.B4.NIR = init(rL, fun=runif)
rL.B3.Red = init(rL, fun=runif)          

# NDVI for 2007
# .............
rL.B4.NIR = wofTara.2007.LSSR.rB[["Band4.NIR"]]
rL.B3.Red = wofTara.2007.LSSR.rB[["Band3.Red"]]
wofTara.2007.NDVI.rL.3b2 = overlay(rL.B4.NIR, rL.B3.Red, fun=calc.NDVI.f2)
wofTara.2007.NDVI.rL.3b2

## class       : RasterLayer 
## dimensions  : 408, 451, 184008  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -0.08312343, 0.7602862  (min, max)

# NDVI for 2017
# .............
rL.B4.NIR = wofTara.2017.LSSR.rB[["Band4.NIR"]]
rL.B3.Red = wofTara.2017.LSSR.rB[["Band3.Red"]]
wofTara.2017.NDVI.rL.3b2 = overlay(rL.B4.NIR, rL.B3.Red, fun=calc.NDVI.f2)
wofTara.2017.NDVI.rL.3b2

## class       : RasterLayer 
## dimensions  : 408, 451, 184008  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -0.6510903, 0.8544872  (min, max)

#=========================================================================================
# Compare Results (those we can compare)
#=========================================================================================

# NDVI for 2007
# .............
all.equal(wofTara.2007.NDVI.rL.1, wofTara.2007.NDVI.rL.2)

## [1] TRUE

all.equal(wofTara.2007.NDVI.rL.1, wofTara.2007.NDVI.rL.3a2)

## [1] TRUE

all.equal(wofTara.2007.NDVI.rL.1, wofTara.2007.NDVI.rL.3b2)

## [1] TRUE

# NDVI for 2017
# .............
all.equal(wofTara.2017.NDVI.rL.1, wofTara.2017.NDVI.rL.2)

## [1] TRUE

all.equal(wofTara.2017.NDVI.rL.1, wofTara.2017.NDVI.rL.3a2)

## [1] TRUE

all.equal(wofTara.2017.NDVI.rL.1, wofTara.2017.NDVI.rL.3b2)

## [1] TRUE

# Remove unneeded Raster Objects
# ------------------------------
ls(pattern="rL.2|rL.3")

## [1] "wofTara.2007.NDVI.rL.2"   "wofTara.2007.NDVI.rL.3a2"
## [3] "wofTara.2007.NDVI.rL.3b2" "wofTara.2017.NDVI.rL.2"  
## [5] "wofTara.2017.NDVI.rL.3a2" "wofTara.2017.NDVI.rL.3b2"

rm(list=ls(pattern="rL.2|rL.3"))

We now plot the NDVI rasters using the levelplot function from the rasterVis package, which produces nice plots with more information than the standard plot function in the raster package`. We plot low NDVI values in red, intermediate values in yellow, and high values in green (the higher NDVI the darker the green).

To plot the NDVI’s of both year we first create a RasterBrick with both RasterLayers. We could plot both RasterLayers side to side but by default this would produce smaller plots as the plots would include density plots on their margins, which are not strictly comparable. They would also include two scales; plotting a RasterBrick a single scale is shown.

#=========================================================================================
# Plot NDVI rasters
#=========================================================================================

# Create the color palette
# ------------------------ 
# Red: low NDVIs, Yellow: middle NDVIs; Green: high NDVIs (Darker Greens indicate Higher NDVIs)
NDVI.breaks = seq(0,1,by=0.01)
NDVI.cols = colorRampPalette(c("red", "yellow", "darkgreen"))(length(NDVI.breaks)-1)

# Create side by side plots
# -------------------------
wofTara.NDVIs.rB = brick(wofTara.2007.NDVI.rL.1, wofTara.2017.NDVI.rL.1)
names(wofTara.NDVIs.rB) = c("NDVI.Autumn 2007", "NDVI.Autumn 2017")
levelplot(wofTara.NDVIs.rB, at=NDVI.breaks, col.regions=NDVI.cols, main="NDVI", names.attr=c("Autumn 2007","Autumn 2017"))

NDVI Change Analysis

Now we explore the change in NDVI in the 10 year period (2007 to 2017). First, we compute the change in NDVI and then we visualise it.

NDVI Change Calculations

To explore the change in NDVI we substract the NDVI values of 2007 from 2017. Therefore, positive values indicate an increase in NDVI, 0 no change, and negative values a decrease in NDVI.

Differences in NDVI can be due to many factors. Typically, the largest difference is due to climatic effects. Here the overall larger NDVIs in 2017 than in 2007 are due to a much wetter start of the year in 2017 than in in 2007. However, there are also specific areas of change. There are more watering points (indicated by very low NDVI values) and there is more coal seam gas infrastructure.

To make changes in NDVI slightly easier to interpret we normalise NDVI values (i.e. (NDVI - mean(NDVI)) / std(NDVI) ), so that 0 represent the mean NDVI and changes are in standard deviations of NDVI. To compute the summary statics values of a raster object (e.g. mean and standard deviation in our case) we need to use the cellStats function in the raster package. See TERN’s DSDP Tutorial ‘Using Raster Data in R’ for details.

As when we calculated NDVI, NDVI changes (both in absolute NDIV values and in stdev. units) can be computed using:

‘Raster Algebra’ (= ‘Raster Maths’): either directly using a formula, or placing this formula in a function and calling the function.
‘Higher Level Functions’, such as calc and overlay.

For normalising the NDVI values in a raster we can use the function normImage in the package RStoolbox. We could incorporate this function both in ‘Raster Algebra’ formula/function and in the function called by the ‘Higher Level Function’. Here we will use ‘Raster Algebra’ both specifying the complete formula and taking advantage of the function normImage.

# Difference in 'raw' NDVI values
# ===============================
wofTara.NDVI.Diff.2017to2007.rL = wofTara.2017.NDVI.rL.1 - wofTara.2007.NDVI.rL.1
wofTara.NDVI.Diff.2017to2007.rL

## class       : RasterLayer 
## dimensions  : 408, 451, 184008  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -1.229459, 0.606049  (min, max)

# Difference in normalise NDVI values
# ===================================

# Method 1: Specifiying the whole formula
# ---------------------------------------
wofTara.NDVI.NormDiff.2017to2007.rL.m1 = 
    ((wofTara.2017.NDVI.rL.1 - cellStats(wofTara.2017.NDVI.rL.1,mean)) / cellStats(wofTara.2017.NDVI.rL.1,sd)) - 
    ((wofTara.2007.NDVI.rL.1 - cellStats(wofTara.2007.NDVI.rL.1,mean)) / cellStats(wofTara.2007.NDVI.rL.1,sd))  
wofTara.NDVI.NormDiff.2017to2007.rL.m1

## class       : RasterLayer 
## dimensions  : 408, 451, 184008  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -10.9933, 3.723755  (min, max)

# Method 2: Using the 'normImage' function in the 'RStoolbox' package
# -------------------------------------------------------------------
wofTara.NDVI.NormDiff.2017to2007.rL.m2 = normImage(wofTara.2017.NDVI.rL.1) - 
                                         normImage(wofTara.2007.NDVI.rL.1)
wofTara.NDVI.NormDiff.2017to2007.rL.m2

## class       : RasterLayer 
## dimensions  : 408, 451, 184008  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1793835, 1807365, -3066365, -3054125  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -10.9933, 3.723755  (min, max)

# Compare results by both Normalisation Methods
# ---------------------------------------------
all.equal(wofTara.NDVI.NormDiff.2017to2007.rL.m1, wofTara.NDVI.NormDiff.2017to2007.rL.m2)

## [1] TRUE

NDVI Change Visualizations

To visualise the changes in ‘raw’ and normalised NDVI we will examine:

The NDVI values of individual years using:
Boxplots: Using the function bwplot from the package rasterVis. This function creates violin plots, a boxplot with a kernel density plot.
Matrix of scatterplots with a fitted loess to investigate the relationship between both years.
The NDVI values of individual years and their change using Density plots: We will build the density plots in two ways:
- Using the function densityplot from the package rasterVis: Simple easy plots.
- Using the function ggplot from the package ggplot2: This looks prettier and is more flexible, but requires more work.
The NDVI change between years by Plotting a RasterLayer of this change (in both scales: raw NDVIs and normalised NDVIs).

We will start by creating a RasterBrick object to pass to many of the plotting functions (e.g. the functions in the package rasterVis).

# =========================
# Create RasterBrick object
# =========================

# 'Raw' NDVIs: Already done above (see 'Plot NDVI rasters' section)
# ...........
wofTara.RawNDVIs.rB = brick(wofTara.2007.NDVI.rL.1, wofTara.2017.NDVI.rL.1) 
names(wofTara.RawNDVIs.rB) = c("Raw.NDVI.2007", "Raw.NDVI.2017")

# Normalised NDVIs
# ................
wofTara.2007.NormNDVI.rL.1 = normImage(wofTara.2007.NDVI.rL.1)
wofTara.2017.NormNDVI.rL.1 = normImage(wofTara.2017.NDVI.rL.1)
wofTara.NormNDVIs.rB = brick(wofTara.2007.NormNDVI.rL.1, wofTara.2017.NormNDVI.rL.1)
names(wofTara.NormNDVIs.rB) = c("Norm.NDVI.2007", "Norm.NDVI.2017")


# ========
# Boxplots
# ========

# Change in 'Raw' NDVI
# ....................
RawNDVI.IndvYrs.bwplot = bwplot(wofTara.RawNDVIs.rB, main="Raw NDVIs")

# Change in Normalised NDVI
# .........................
NormNDVI.IndvYrs.bwplot = bwplot(wofTara.NormNDVIs.rB, main="Normalised NDVIs")


# ========================
# Scatter plots with loess
# ========================

# Change in 'Raw' NDVI
# ....................
RawNDVI.IndvYrs.splot = splom(wofTara.RawNDVIs.rB, main="Raw NDVIs", plot.loess=TRUE, xlab='')

# Change in Normalised NDVI
# .........................
NormNDVI.IndvYrs.splot = splom(wofTara.NormNDVIs.rB, main="Normalised NDVIs", plot.loess=TRUE, xlab='')


# -------------------------------------
# Combine the 4 Plots in a single Graph
# -------------------------------------
grid.arrange( RawNDVI.IndvYrs.bwplot, NormNDVI.IndvYrs.bwplot, 
              RawNDVI.IndvYrs.splot, NormNDVI.IndvYrs.splot,
              nrow=2, top="NDVI - Individuals Years" )

# =============    
# Density Plots
# =============

# ----------------------------------------------------------
# Method 1: Using 'densityplot' from the package 'rasterVis`
# ----------------------------------------------------------

# Change in 'Raw' NDVI
# ....................
RawNDVI.IndvYrs.dp = densityplot(wofTara.RawNDVIs.rB, main="Individual Raw NDVIs")
RawNDVI.Diff.dp = densityplot(wofTara.NDVI.Diff.2017to2007.rL, main="Change in Raw NDVIs (2007-2017)")

# Change in Normalised NDVI
# .........................
NormNDVI.IndvYrs.dp = densityplot(wofTara.NormNDVIs.rB, main="Individual Normalised NDVIs")
NormNDVI.Diff.dp = densityplot(wofTara.NDVI.NormDiff.2017to2007.rL.m2, main="Change in Normalised NDVIs (2007-2017)")

# Combine 4 Plots in a Graph
# ..........................
grid.arrange( RawNDVI.IndvYrs.dp, RawNDVI.Diff.dp, NormNDVI.IndvYrs.dp, NormNDVI.Diff.dp,
              nrow=2, top="NDVI: Individuals Years & Change" )

# ---------------------------------------------------
# Method 2: Using 'ggplot' from the package 'ggplot2`
# ---------------------------------------------------

# Create a DataFrame with the required values 
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Create a DataFrame with the required values to pass to function 'ggplot' in the package 'ggplot2' 
# Required values: NDVIs for 2007, 2017, and their difference in both scales: raw and normalised.
RawNDVI.2007 = values(wofTara.2007.NDVI.rL.1)
RawNDVI.2017 = values(wofTara.2017.NDVI.rL.1)
RawNDVI.Diff = values(wofTara.NDVI.Diff.2017to2007.rL)
NormNDVI.2007 = values(wofTara.2007.NormNDVI.rL.1)
NormNDVI.2017 = values(wofTara.2017.NormNDVI.rL.1)
NormNDVI.Diff = values(wofTara.NDVI.NormDiff.2017to2007.rL.m2)
NDVIs.df = data.frame( RawNDVI.2007, RawNDVI.2017, RawNDVI.Diff, 
                       NormNDVI.2007, NormNDVI.2017, NormNDVI.Diff )
#summary(NDVIs.df)
                               

# Change in 'Raw' NDVI
# ~~~~~~~~~~~~~~~~~~~~

# Density plot of Raw NDVI values for Individual Years: 
# .....................................................
# Subset Dataset and convert wide format to long format
RawNDVI.IndvYrs.dflf = melt(NDVIs.df[c("RawNDVI.2007","RawNDVI.2017")])
RawNDVI.IndvYrs.dp2 = ggplot(RawNDVI.IndvYrs.dflf, aes(x=value, fill=variable)) +
                      geom_density(alpha=0.25)
                   
# Density plot of the change in Raw NDVI values
# .............................................
RawNDVI.Diff.dp2 = ggplot(NDVIs.df, aes(x=RawNDVI.Diff)) + 
                   geom_density(alpha=0.25,color="darkgreen", fill="lightgreen")
                   

# Change in Normalised NDVI
# ~~~~~~~~~~~~~~~~~~~~~~~~~

# Density plot of Raw NDVI values for Individual Years: 
# .....................................................
# Subset Dataset and convert wide format to long format
NormNDVI.IndvYrs.dflf = melt(NDVIs.df[c("NormNDVI.2007","NormNDVI.2017")])
NormNDVI.IndvYrs.dp2 = ggplot(NormNDVI.IndvYrs.dflf, aes(x=value, fill=variable)) +
                       geom_density(alpha=0.25)
                   
# Density plot of the change in Raw NDVI values
# .............................................
NormNDVI.Diff.dp2 = ggplot(NDVIs.df, aes(x=NormNDVI.Diff)) + 
                    geom_density(alpha=0.25,color="darkgreen", fill="lightgreen")

                    
# Combine the 4 Plots in a single Graph
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
grid.arrange( RawNDVI.IndvYrs.dp2, RawNDVI.Diff.dp2, NormNDVI.IndvYrs.dp2,   NormNDVI.Diff.dp2,
              nrow=2, top="NDVI: Individuals Years & Change" )

# ============================================================
# Plot of the Raster Layer with the Normalised Change in NDVI
# ============================================================
# They have very differnt scales (raw NDVI and normalised NDVI values)

# -------------------
# Change in Raw NDVI 
# -------------------
#wofTara.NDVI.Diff.2017to2007.rL

# Create the color palette
# ........................
abs.max = max(abs(values(wofTara.NDVI.Diff.2017to2007.rL)), na.rm=TRUE)
# Red: low NDVIs, Yellow: middle NDVIs; Green: high NDVIs (Darker Greens indicate Higher NDVIs)
RawNDVI.Diff.breaks = seq(-abs.max,abs.max, by=0.01)
RawNDVI.Diff.cols = colorRampPalette(c("red", "yellow", "darkgreen"))(length(RawNDVI.Diff.breaks)-1)

# Create Raster Plots with Density Plots
# ......................................
wofTara.RawNDVIDiff.p = levelplot( wofTara.NDVI.Diff.2017to2007.rL, at=RawNDVI.Diff.breaks, 
                                   col.regions=RawNDVI.Diff.cols, main="Raw NDVI Change (2017 - 2007)" )                               
                                   
# -------------------------                                
# Change in Normalised NDVI
# -------------------------
#wofTara.NDVI.NormDiff.2017to2007.rL.m2

# Create the color palette
# ........................
abs.max = max(abs(values(wofTara.NDVI.NormDiff.2017to2007.rL.m2)), na.rm=TRUE)
# Red: low NDVIs, Yellow: middle NDVIs; Green: high NDVIs (Darker Greens indicate Higher NDVIs)
NormNDVI.Diff.breaks = seq(-abs.max,abs.max, by=0.01)
NormNDVI.Diff.cols = colorRampPalette(c("red", "yellow", "darkgreen"))(length(NormNDVI.Diff.breaks)-1)

                            
# Create Raster Plots with Density Plots
# ......................................
wofTara.NormNDVIDiff.p = levelplot( wofTara.NDVI.NormDiff.2017to2007.rL.m2,   at=NormNDVI.Diff.breaks, col.regions=NormNDVI.Diff.cols, main="Normalised NDVI Change (2017 - 2007)")
                            
# -----------------                         
# Plot Raster Plots
# -----------------     
# `grid.arrange` provides multiple plots per page. They look good on screen, 
# but might be too small in the Workshop/Tutorial report. 
grid.arrange(wofTara.RawNDVIDiff.p, wofTara.NormNDVIDiff.p, nrow=1) # Too Small for Tut.

#wofTara.RawNDVIDiff.p   # Used if 'grid.arrange' leads to plots too small
#wofTara.NormNDVIDiff.p  # Used if 'grid.arrange' leads to plots too small

Examine what’s happening in the Area of Most Noticeable Greening up

Aside from areas of where infrastructure has been put in and a small patch of clearing in the south-east corner, there is an area of noticeable greening up in the north-east.

To examine what is happening here overtime we can plot a time series of the green fraction. The green fraction is one of the bands of the fractional cover product. This is another AusCover product hosted by TERN.

Then we fit a trend line to the green data over the period of interest.

The server hosts time series stacks of cover fractions in virtual files that make pulling out time series a breeze. However, individual bands could also be iterated without too many problems.

Greening Area: Define its Extent & Find out where it is

To find our bearings we start by displaying the Greening Area. We do this in different ways: * Download and Display a map of the Greening Area from Stamen Maps. * Plot Satellite images of our study area (west of Tara) for years 2007 & 2017 with a (red) box indicating the Greening Area. * Plot the Normalised NDVI Change image with a (blue) box indicating the Greening Area.

The plots for the Satellite (Surface Reflectance) images are ggplot2 plots. The Normalised NDVI Change images is a lattice plot. The procedure to add a polygon box to ggplot2 & lattice plots differ. We plot both types of plots to show how to add a polygon in each case.

# Define the Extent of the Area to Explore
# ========================================
GreeningArea.extent = extent(150.4436, 150.4458, -27.0062, -27.0037)
GreeningArea.extent

## class       : Extent 
## xmin        : 150.4436 
## xmax        : 150.4458 
## ymin        : -27.0062 
## ymax        : -27.0037

# Plot Web Map for the Greening Area
# ==================================

# Get Map 
# -------
GreeningArea.StamenMap = get_stamenmap(GreeningArea.extent[c(1,3,2,4)], maptype="terrain")

# Create Plot of the Map
# ----------------------
# Stamen map of this area is not very detailed
GreeningArea.P1 = 
ggmap(GreeningArea.StamenMap) + labs(title= "Greening Area", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"), 
      axis.title = element_text(size=11, face="bold"), axis.text=element_text(size=9) )
GreeningArea.P1

# Plot Greening Area over Satellite and NDVI Change Images
# ========================================================

# Create Polygon with the Greening Area Extent
# --------------------------------------------
GreeningArea.SP = as(GreeningArea.extent, "SpatialPolygons")
GreeningArea.SP  # Has no CRS

## class       : SpatialPolygons 
## features    : 1 
## extent      : 150.4436, 150.4458, -27.0062, -27.0037  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

proj4string(GreeningArea.SP) = CRS("+init=epsg:4326") # Add CRS
GreeningArea.SP

## class       : SpatialPolygons 
## features    : 1 
## extent      : 150.4436, 150.4458, -27.0062, -27.0037  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

GreeningArea.EPSG3577.SP = spTransform(GreeningArea.SP, crs("+init=epsg:3577"))
GreeningArea.EPSG3577.SP

## class       : SpatialPolygons 
## features    : 1 
## extent      : 1801860, 1802113, -3057837, -3057529  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:3577 +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

# Plot Greening Area on Satellite Images
# --------------------------------------
# Plots 'wofTara.2007.LSSR.p' and 'wofTara.2017.LSSR.p' were created above.

# 2007 Image
# ~~~~~~~~~~
GreeningArea.2007.p = wofTara.2007.LSSR.p + 
                      geom_polygon(data=GreeningArea.EPSG3577.SP, aes(x=long, y=lat), 
                                   color="red", alpha=0)

# 2007 Image
# ~~~~~~~~~~
GreeningArea.2017.p = wofTara.2017.LSSR.p + 
                      geom_polygon(data=GreeningArea.EPSG3577.SP, aes(x=long, y=lat), 
                                   color="red", alpha=0)

# Plot Satellite Images from 2007 & 2017 with a Red Box for the Greening Area
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# `grid.arrange` provides multiple plots per page. They look good on screen, 
# but might be too small in the Workshop/Tutorial report. 
grid.arrange(GreeningArea.2007.p, GreeningArea.2017.p, nrow=1)

#GreeningArea.2007.p   # Used if 'grid.arrange' leads to plots too small
#GreeningArea.2017.p   # Used if 'grid.arrange' leads to plots too small

# Plot Greening Area on NDVI Change Image
# ---------------------------------------
# Parameters 'NormNDVI.Diff.breaks' & 'NormNDVI.Diff.cols' were created above 
# (when genrating plot 'wofTara.NormNDVIDiff.p'.

levelplot( wofTara.NDVI.NormDiff.2017to2007.rL.m2, at=NormNDVI.Diff.breaks, 
           col.regions=NormNDVI.Diff.cols, 
           main="Normalised NDVI Change (2017 - 2007) with Greening Area in blue box") +
layer(sp.polygons(GreeningArea.EPSG3577.SP, col="blue", fill=NA))

Download the dataset with the Green Cover Fraction

We start by downloading the vrt file containing the green cover fraction. VRT is format (Virtual Format) driver for GDAL (Geospatial Data Abstraction Library) that allows building virtual GDAL datasets from other GDAL datasets. ‘.vrt’ files contain ‘VRT’ descriptions of datasets typically saved in XML format. For further details see ‘GDAL Virtual Format Tutorial’ by GDAL organization here

NOTE: That the name of the vrt file might have changed since this document was written, as files get updated (i.e. they grow) with new data.

# Data Path and Name
data.path = "http://qld.auscover.org.au/public/data/landsat/seasonal_fractional_cover/fractional_cover/aus"
#data.fn = "lztmre_aus_s198712201705_dima2_green.vrt"
#data.fn = "lztmre_aus_s198712201805_dima2_green.vrt"
data.fn = "lztmre_aus_s198712201808_dima2_green.vrt" # It's changed again!!!

# Download the data file (if it doesn't work try method='wget' or functions in Library 'RCurl')
download.file(url=paste(data.path, data.fn, sep="/"), destfile='GreenCoverFraction.vrt', method='auto')
list.files()

##  [1] "ESA18_Workshop_AusCover_Overview.docx"   
##  [2] "ESA18_Workshop_AusCover_P1.html"         
##  [3] "ESA18_Workshop_AusCover_P1.R"            
##  [4] "ESA18_Workshop_AusCover_P1.Rmd"          
##  [5] "ESA18_Workshop_AusCover_P1_NoStamen.html"
##  [6] "ESA18_Workshop_AusCover_P1_NoStamen.R"   
##  [7] "ESA18_Workshop_AusCover_P1_NoStamen.Rmd" 
##  [8] "ESA18_Workshop_AusCover_P2.html"         
##  [9] "ESA18_Workshop_AusCover_P2.R"            
## [10] "ESA18_Workshop_AusCover_P2.Rmd"          
## [11] "ESA18_Workshop_AusCover_P2_files"        
## [12] "ESA18_Workshop_AusCover_P2_NoStamen.html"
## [13] "ESA18_Workshop_AusCover_P2_NoStamen.R"   
## [14] "ESA18_Workshop_AusCover_P2_NoStamen.Rmd" 
## [15] "ESA18_Workshop_AusCover_P2_v6f.html"     
## [16] "ESA18_Workshop_AusCover_P2_v6f.R"        
## [17] "ESA18_Workshop_AusCover_P2_v6f.Rmd"      
## [18] "ESA18_Workshop_AusCover_ReadMe.docx"     
## [19] "GreenCoverFraction.vrt"                  
## [20] "LSSR_2007.vrt"                           
## [21] "LSSR_2017.vrt"

Extract and Prepare (Crop and Re-code) the Bands of Interest (2007 to 2017)

The green cover fraction vrt file currently contains 122 bands/layers. Four new bands are created and added each year. Bands 77-120 correspond to the period of interest, 2007-2017 (one band has already been added for 2018). You can check this (and the file structure) opening it in text editor such as Notepad.

We create a new function to Load, Subset (i.e. Crop), Recode (the no-data values) and Combine the bands of interest. This function is partially similar to the function ‘get.B3B4B5.rB.f’ that we created above. Then we call this function to create a RasterBrick containing the Green Cover Fraction bands of interest (bands 77-120 for the period 2007-2017) subsetted to the extent of the Greening Area we want to explore.

#-----------------------------------------------------------------------------------------
# Create a function to Extract, Subset, Recode, and Combine the Bands into a Raster Brick
#-----------------------------------------------------------------------------------------

getprep.Bands.rB.f = function(dl.fn, bands.v, rB.extent, rB.crs){

    # Create cookie cutter: Extent projected to CRS=EPSG:3577
    # =======================================================
    #print(rB.extent)
    # Create a Polygon from Extent, Reproject the Poloygon and use its extension
    rB.extent.GeoCoord.SP = as(rB.extent, "SpatialPolygons")
    #print(rB.extent.GeoCoord.SP)  # Has no CRS
    proj4string(rB.extent.GeoCoord.SP) = CRS("+init=epsg:4326")
    #print(rB.extent.GeoCoord.SP)
    rB.extent.reprj.SP = spTransform(rB.extent.GeoCoord.SP, rB.crs)
    #print(rB.extent.reprj.SP)
    rB.extent.reprj = extent(rB.extent.reprj.SP)
    #print(rB.extent.reprj)

    # For each Band
    # =============
    for ( band.cnt in bands.v ) {
    
        #print(band.cnt)
    
        # Load the data to RasterLayer object
        # -----------------------------------
        GFC.rL = raster(dl.fn, band=band.cnt)
        #print(GFC.rL)
    
        # Subset (i.e. crop) raster layer to area of interest
        # ---------------------------------------------------
        GFC.Subset.rL = crop(GFC.rL, rB.extent.reprj)
        is.na(GFC.Subset.rL)
        #print(GFC.Subset.rL)
    
        # Replace values < 0 with NAs
        # ---------------------------
        # GFC.Subset.rL[GFC.Subset.rL < 0] = NA
    
        # Add Raster Layer to Raster Brick
        # --------------------------------
        if ( band.cnt == bands.v[1] )
            GFC.Subset.rB = brick(GFC.Subset.rL)
        else
            GFC.Subset.rB = addLayer(GFC.Subset.rB, GFC.Subset.rL)
            
    } # for ( band.cnt in bands.v ) {

    # Re-conver to RasterBrick object (if prefered)
    # =============================================
    # 'addLayer' (used in the for loop) returns a RasterStack object, so we need to re-convert 
    # object to a RasterBrick (if this is the prefered object class).
    GFC.Subset.rB = brick(GFC.Subset.rB)
    
    # Add Layer Names to Raster Brick
    # ===============================
    
    # Create a vector of Raster Layer names (including the Band, Year, and Quarter)
    band.part = paste("Band", bands.v, sep="")
    initial.year = 1988
    date.part = as.character(formatC(initial.year+bands.v/4-0.25,digits=2,format="f")) # '-0.25 'cos statrs in band1 not band0
    band.names = paste(band.part, date.part, sep="_")
    
    # Add the band names to the RasterBrick
    names(GFC.Subset.rB) = band.names
    #print(GFC.Subset.rB)

    # Return RasterBrick with the Required Bands Subset to the desired Extent
    # =======================================================================
    return(GFC.Subset.rB)   
 
} # getprep.Bands.rB.f = function(dl.fn, bands.v, rB.extent, rB.crs){


#----------------------------------------------------------------------------------------
# Call the function to Create a RasterBrick with processed Bands for period 2007 to 2017
#----------------------------------------------------------------------------------------

GreeningArea.rB =getprep.Bands.rB.f( dl.fn='GreenCoverFraction.vrt', bands.v=77:120, 
                                     rB.extent=GreeningArea.extent, rB.crs=CRS("+init=epsg:3577") )
GreeningArea.rB

## class       : RasterBrick 
## dimensions  : 11, 9, 99, 44  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 1801845, 1802115, -3057845, -3057515  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : Band77_2007.00, Band78_2007.25, Band79_2007.50, Band80_2007.75, Band81_2008.00, Band82_2008.25, Band83_2008.50, Band84_2008.75, Band85_2009.00, Band86_2009.25, Band87_2009.50, Band88_2009.75, Band89_2010.00, Band90_2010.25, Band91_2010.50, ... 
## min values  :            100,            100,            100,            116,            121,            112,            109,            112,            113,            116,            120,            110,            118,            126,            111, ... 
## max values  :            124,            128,            138,            149,            150,            142,            136,            128,            127,            135,            145,            127,            140,            153,            146, ...

Compute the Mean Green Cover Fraction in the Greening Area for each raster layer

First we compute a vector with the mean Green Cover Fraction (GCF) in the Greening Area (GA) for each raster layers. The raster layer values are in scale 0-255, so we re-scale the mean value to scale of 0-100 to represent the Green Cover Fraction in Percentage.

Then we create a vector with the Times (Years and Quarters) for each of these raster layers. The function that we used to create the Raster Brick with all the relevant data computed the date for each layer and assigned it as part of the corresponding band (i.e. layer) name. Here, we take advantage of this and extract the Times from the bands names. We could also have calculated it in a similar fashion as it was done within the function.

# Compute the Mean GCF in the Greening Area for each of the Bands
# ================================================================
GCF.GA.Means = cellStats(GreeningArea.rB*100/255, mean)
GCF.GA.Means

##  Band77_2007.00  Band78_2007.25  Band79_2007.50  Band80_2007.75 
##        42.07170        43.42246        45.26837        52.84611 
##  Band81_2008.00  Band82_2008.25  Band83_2008.50  Band84_2008.75 
##        53.65815        48.69479        47.51040        46.94791 
##  Band85_2009.00  Band86_2009.25  Band87_2009.50  Band88_2009.75 
##        46.27055        48.67895        52.75500        45.90216 
##  Band89_2010.00  Band90_2010.25  Band91_2010.50  Band92_2010.75 
##        49.40780        55.65855        50.83779             NaN 
##  Band93_2011.00  Band94_2011.25  Band95_2011.50  Band96_2011.75 
##        59.07704        60.67340        52.73123        52.49752 
##  Band97_2012.00  Band98_2012.25  Band99_2012.50 Band100_2012.75 
##             NaN        55.62370        54.23648        49.41573 
## Band101_2013.00 Band102_2013.25 Band103_2013.50 Band104_2013.75 
##             NaN        57.22321        52.92533        50.03763 
## Band105_2014.00 Band106_2014.25 Band107_2014.50 Band108_2014.75 
##        51.25767        56.03882        54.03446        49.59794 
## Band109_2015.00 Band110_2015.25 Band111_2015.50 Band112_2015.75 
##        53.66211        54.51773        54.23252        51.46762 
## Band113_2016.00 Band114_2016.25 Band115_2016.50 Band116_2016.75 
##        54.96534        57.83720        56.28045        54.83462 
## Band117_2017.00 Band118_2017.25 Band119_2017.50 Band120_2017.75 
##        55.67835        58.90275        56.31214        52.38661

# Create a vector with the Times for each of Bands 
# ================================================
#names(GCF.GA.Means)
GCF.GA.Times = as.numeric( sub(".*_","", names(GCF.GA.Means)) )
# OR:  as.numeric( sapply(strsplit(names(GCF.GA.Means),split="[_]"),"[",2) )
GCF.GA.Times

##  [1] 2007.00 2007.25 2007.50 2007.75 2008.00 2008.25 2008.50 2008.75
##  [9] 2009.00 2009.25 2009.50 2009.75 2010.00 2010.25 2010.50 2010.75
## [17] 2011.00 2011.25 2011.50 2011.75 2012.00 2012.25 2012.50 2012.75
## [25] 2013.00 2013.25 2013.50 2013.75 2014.00 2014.25 2014.50 2014.75
## [33] 2015.00 2015.25 2015.50 2015.75 2016.00 2016.25 2016.50 2016.75
## [41] 2017.00 2017.25 2017.50 2017.75

Plot the Change in Green Cover Fraction in the Greening Area from 2007 to 2017

To explore the Change in Green Cover Fraction in the Greening Area in the study period (2007-2017) we create two types of plots:

Scatterplot with a trend: To show the trend we use different 3 types of smooth: (1) a linear model fit, (2) a linear model fit with 2 degree polynomial (using the function poly), and (3) a locally weighted regression (loess). For prettier and more flexible graphs we use the ggplot function from the ggplot2` package.
Violin Plots: Violin plots combine Boxplots and Density plots (see above). We use the function bwplot from the package rasterVis to create this type of plot.

# Scatter-plot (with lines) using different Smooths
# =================================================
# 'ggplot' requires a dataframe, so we first create a dataframe containing the 
# Green Cover Fractions (%) and Times (Years & Quarters).

# Create dataframe
# ----------------
GCF.GA.df = data.frame(Time=GCF.GA.Times, Mean=GCF.GA.Means)
summary(GCF.GA.df)

##       Time           Mean      
##  Min.   :2007   Min.   :42.07  
##  1st Qu.:2010   1st Qu.:49.42  
##  Median :2012   Median :52.85  
##  Mean   :2012   Mean   :52.35  
##  3rd Qu.:2015   3rd Qu.:55.62  
##  Max.   :2018   Max.   :60.67  
##                 NA's   :3

# Create base plot (base + points)
# --------------------------------
base.p = ggplot(GCF.GA.df, aes(x=Time, y=Mean)) + 
         theme(plot.title=element_text(hjust = 0.5, size=14, face="bold")) +
         labs(x="Time (Year.Quarter)", y="Green Cover Fraction (%)")    + 
         scale_x_continuous(breaks=2007:2018) +
         geom_line(color="blue")
#base.p

# Add different types of Smooths to the Base Plot
# -----------------------------------------------

# No Trend
base.NoSmooth.p = base.p + ggtitle("Green Cover Fraction in Greening Area (2007 to 2017)")
#base.NoSmooth.p

# Linear Model (LM) Fit
base.SmoothLM.p = base.p + ggtitle("Green Cover Fraction + LM Trend") + 
                  stat_smooth(method = "lm", formula = y ~ x, size = 1)
#base.SmoothLM.p

# LM with 2nd Degree Polynomial Fit (using the function 'poly')
base.SmoothPoly2D.p = base.p + 
                      ggtitle("Green Cover Fraction + LM with 2 Deg. Poly. Trend") + 
                      stat_smooth(method = "lm", formula = y ~ poly(x,2), size = 1)
#base.SmoothPoly2D.p

# Logically Weighted Regression (Loess)
base.SmoothLoess.p = base.p + ggtitle("Green Cover Fraction + LOESS Trend") + 
                     stat_smooth(method = "loess", formula = y ~ x, size = 1)
#base.SmoothLoess.p

# Display all Scatter-plots (with Lines) with different types of Smooths
# ----------------------------------------------------------------------
# `grid.arrange` provides multiple plots per page. They look good on screen, 
# but might be too small in the Workshop/Tutorial report. 
grid.arrange(base.NoSmooth.p, base.SmoothLM.p, base.SmoothPoly2D.p, base.SmoothLoess.p,
                  nrow=2)

#base.NoSmooth.p        # Used if 'grid.arrange' leads to plots too small
#base.SmoothLM.p        # Used if 'grid.arrange' leads to plots too small
#base.SmoothPoly2D.p    # Used if 'grid.arrange' leads to plots too small
#base.SmoothLoess.p # Used if 'grid.arrange' leads to plots too small


# Violin Plots (Boxplot + Density Plots)
# ======================================
# Rather than using the default tick labels (too long including the Band and Time), we 
# use tick labels including only the Time (Year & Quarter).

bwplot( GreeningArea.rB*100/255, 
        main="Green Cover Fraction in the Greening Area (2007 to 2017)", 
            xlab="Time (Year.Quarter)", ylab="Green Cover Fraction (%)", 
            scales=list(x=list(labels=format(GCF.GA.Times,nsmall=2),rot=45)) )

ESA18 TERN Workshop: Landscape Assessment (AusCover) - Remote Sensing - (Part 2)

"Simple Change Detection in Raster using TERN web services"

NORMALIZED DIFFERENCE VEGETATION INDEX (NDVI) CHANGE ANALYSIS

Preparation: Getting ready

Explore an Area where NDVI Change might have occurred

Create a function to Download, Read, Subset, and Combine the Landsat Surface Reflectance Data.

Obtain Landsat Surface Reflectance (LSSR) Data for the required Extent and Period using the new function (`get.B3B4B5.rB.f`) & Plot them

Compute NDVI from the extracted Subsets

NDVI Change Analysis

NDVI Change Calculations

NDVI Change Visualizations

Examine what’s happening in the Area of Most Noticeable Greening up

Greening Area: Define its Extent & Find out where it is

Download the dataset with the Green Cover Fraction

Extract and Prepare (Crop and Re-code) the Bands of Interest (2007 to 2017)

Compute the Mean Green Cover Fraction in the Greening Area for each raster layer

Plot the Change in Green Cover Fraction in the Greening Area from 2007 to 2017

ESA18 TERN Workshop: Landscape Assessment (AusCover) - Remote Sensing - (Part 2)

"Simple Change Detection in Raster using TERN web services"

NORMALIZED DIFFERENCE VEGETATION INDEX (NDVI) CHANGE ANALYSIS

Preparation: Getting ready

Explore an Area where NDVI Change might have occurred

Create a function to Download, Read, Subset, and Combine the Landsat Surface Reflectance Data.

Obtain Landsat Surface Reflectance (LSSR) Data for the required Extent and Period using the new function (get.B3B4B5.rB.f) & Plot them

Compute NDVI from the extracted Subsets

NDVI Change Analysis

NDVI Change Calculations

NDVI Change Visualizations

Examine what’s happening in the Area of Most Noticeable Greening up

Greening Area: Define its Extent & Find out where it is

Download the dataset with the Green Cover Fraction

Extract and Prepare (Crop and Re-code) the Bands of Interest (2007 to 2017)

Compute the Mean Green Cover Fraction in the Greening Area for each raster layer

Plot the Change in Green Cover Fraction in the Greening Area from 2007 to 2017

Obtain Landsat Surface Reflectance (LSSR) Data for the required Extent and Period using the new function (`get.B3B4B5.rB.f`) & Plot them